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Ab initio phase diagram of BaTiOs under epitaxial strain revisited 
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We revisit the phase diagram of BaTiC >3 under biaxial strain using a first principles-based effective Hamil¬ 
tonian approach. We show that, in addition to the tetragonal (c), quasi-rhombohedral (r), and quasi- 
orthorhombic ( aa ) ferroelectric phases, that have been discussed previously, there are temperature and strain 
regions, in particular under tensile strain, where the system decomposes into multi-domain structures. In 
such cases, the strained system, at least on a local level, recovers the same phase sequence as the unclamped 
bulk material. Furthermore, we extend these results from the case of “uniform” biaxial strain to the situation 
where the two in-plane lattice constants are strained differently and show that similar considerations apply 
in this case. 


The optimization of ferroelectric materials by epitaxial 
growth and interface-mediated strain is nowadays a well- 
established and highly successful method.— However, the 
experimental determination of strain-temperature phase 
diagrams is quite challenging, since only specific strain 
values, corresponding to the given lattice mismatch with 
a specific substrate, can be investigated. Therefore, the 
theoretical modeling of strain-dependent phase diagrams 
is highly relevant. 

An important case is the prototypical ferroelectric 
BaTiC >3 (BTO), which, in its free bulk form, exhibits 
one paraelectric and three different ferroelectric struc¬ 
tures as function of temperature^ and thus gives rise 
to a rich strain dependence. Different levels of sophis¬ 
tication have been used to model/calculate the strain- 
dependent phase diagram of BTO, however, leading in 
part to conflicting results. First, various calculations 
based on Ginzburg-Landau-Devonshire theory have been 
performed, yielding qualitatively consistent phase dia¬ 
grams as long as only mono-domain phases are taken into 
account—^ Once multi-domain configurations are con¬ 
sidered, different (meta-) stable domain patterns have 
been found, depending on the a priori assumptions of 
the models. 6,7 More recently, phase field simulations have 
been used in order to simulate different phases and do¬ 
main structures without a priori assumptions— 2 

Both approaches, however, require parameters that 
have to be obtained from experimental data, and different 
choices for these parameters can lead to significantly dif¬ 
ferent results— From this perspective, first principles- 
based calculations are very attractive, since in principle 
they do not require fitting. Indeed, first principles-based 
strain-temperature phase diagrams for BTO have been 
calculated^ - — but only mono-domain states have been 
found for simulation cell sizes of about 12-16 unit cells 
along each cartesian direction. 

Apart from determining the stability of different 
phases, the mechanical boundary conditions can also 
modify the domain pattern of the material. For example, 
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the formation of 90° domain walls provides an efficient 
way for tetragonal perovskites to partially relax the elas¬ 
tic energy under clamping to a periodic substrate.— In 
PbTiC> 3 , such domain configurations have been predicted 
theoretically even without epitaxial constraints— In con¬ 
trast, a rather large domain wall energy of 0.9 meV/A 2 
has been obtained for 90° domain walls in BTOA^ and 
it has been predicted that such domains are not sta¬ 
ble for unclamped BTO j 20 i 21 Experimentally, both multi- 
domain and mono-domain phases have been reported for 
BTO grown under tensile strain— — 

In the present paper, we show the presence of multi- 
domain states within a first principles-based effective 
Hamiltonian approach for sufficiently large simulation 
cells. For uniform biaxial tensile strain, the paraelec¬ 
tric phase first transforms to a multi-domain ferroelectric 
phase with local polarization along (100) and under fur¬ 
ther cooling a mono-domain state polarized along (110) 
is obtained. For small strains a multi-domain state with 
both in-plane and out-of-plane polarization is observed, 
while under compressive strain only mono-domain phases 
appear. Qualitatively the same trends are found for non- 
uniform biaxial strain, where a large variety of multi¬ 
domain states exist. 

We use the feram code^i which is based on the effec¬ 
tive Hamiltonian discussed in Refs. Hi], and and 
the parametrization given in Ref. l27l We perform molec¬ 
ular dynamics simulations in the canonical ensemble, em¬ 
ploying a Nose-Poincare thermostat^ periodic boundary 
conditions, and a simulation box of 32x32x32 unit cells, 
i.e. corresponding to periodic images with about 13 nm 
distance. An ideal bulk material is modeled, without 
surfaces and depolarizing fields. To simulate the effect of 
epitaxial strain, the elements 771 and 772 of the homoge¬ 
neous strain tensor (in standard Voigt notation) are fixed 
to the external strain 77, and t)q is set to zero. We define 
77 = 0 corresponding to a lattice constant of 3.996 A, as 
obtained within our approach for the free system directly 
above T q. In the following, the two clamped in-plane di¬ 
rections are denoted as a and b, and the free out-of-plane 
direction as c. 

Fig. [T| shows the calculated temperature dependence of 
the total electric polarization P, separated into in-plane 
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FIG. 1. (Color online) Calculated temperature-dependent po¬ 
larization along (a) [110] and (b) [001] for different values of 
biaxial strain applied in the (001) plane (from —0.85 % to 
+ 1.75 % in steps of 0.1 %). Solid (black) line: unstrained case 
(rj = 0); Dashed (blue) lines: tensile strain; Dash-dotted (red) 
lines: compressive strain. 


and out-of-plane components, for different values of the 
epitaxial strain (ry a = r]b). It can be seen that for zero 
strain both in-plane and out-of-plane components of P 
appear at the same temperature, while under compres¬ 
sive strain the appearance of the out-of-plane (in-plane) 
components are shifted to higher (lower) temperatures, 
and vice versa under tensile strain. This is in agreement 
with Ref. HH and [lg]. However, we also observe that the 
temperature dependence of the in-plane polarization un¬ 
der tensile strain exhibits an unusual “kink-like” feature 
at a (strain-dependent) temperature Tg, somewhat below 
the overall critical temperature Tq- 

Further inspection of the local soft mode configurations 
reveals that in the temperature range between Tc and Tg 
the system exhibits a multi-domain state, with average 
in-plane polarization along [110], but local in-plane po¬ 
larization along [100] and [010], separated by 90° domain 
walls parallel to (110) (see Fig. (2]]a)). Below Tg, the 
system adopts a mono-domain state with both local and 
global in-plane polarization along [110]. The formation 
of a multi-domain state under tensile strain allows the 
material to essentially recover the same phase sequence 
as in the free (unclamped) case, at least on a local level. 
At temperatures immediately below Tc, the tetragonal 
ferroelectric phase with polarization along (100) has the 
lowest free energy in the unclamped case— However, a 
corresponding mono-domain state is strongly disfavored 
by the elastic boundary conditions introduced through 
the epitaxial constraint. Thus, by forming the observed 
multi-domain state, the system can lower its overall free 


energy by achieving (100) polarization locally, at the cost 
of introducing energetically unfavorable domain walls. 
Thereby, the epitaxial constraint r] a = rjb > 0 promotes 
in-plane polarization with an equal volume fraction of 
[100] and [010] domains. However, for cell sizes below 
20x20x20, the energy penalty for forming domain walls 
exceeds the gain in free energy within the (lOO)-polarized 
domains, and consequently the multi-domain phase has 
been overlooked in previous ab initio simulations em¬ 
ploying smaller simulation cells. Under further cooling, 
the (llO)-polarized phase becomes more favorable, cf. 
Ref. |3Q| for the unclamped material, and thus a transi¬ 
tion into a corresponding mono-domain state, which is 
also compatible with the epitaxial constraint under ten¬ 
sile strain, occurs at Tg. 

For local (100) polarization, only 180° domain 
walls and 90° walls parallel to {110} are possible by 
symmetry.— The latter type is indeed observed in our 
simulations (see Fig. [2])a) ), whereas 180° domain walls 
do not relax any elastic energy and are thus not favor¬ 
able in the present case. The polarization across the 90° 
domain wall can be separated into the nearly constant 
polarization perpendicular to the wall ( P ±) and the po¬ 
larization parallel to the wall (P\\), which approximately 
follows a tanh-prohlci 19 i 29 


P\\ (x) = P 0 , || tanh 
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with daw half of the domain wall width, Xq center of one 
wall, To, || > polarization in the domain center, xn width of 
the fit region, and 4 the number of domains in the simula¬ 
tion cell (a minimum of 4 domains is necessary to match 
the domain profile at the periodic boundaries of the sim¬ 
ulation cell). For T = 400K and 0.75% tensile strain, 
the fit shown in Fig.[2](b) yields a domain wall width of 
about 5 nm, in good agreement with literature . 19 ’ 29 

A clear signature of the multi-domain state in the tem¬ 
perature range between Tc and Tg can also be seen in 
various energy contributions, see Fig. [5|c). The mod¬ 
ulation of T|| and the elastic mismatch at the domain 
walls induce an energy penalty in the local mode self 
energy and the inhomogeneous part of the elastic en¬ 
ergy, whereas the coupling energy between local strain 
and soft mode is reduced if the mono-domain state with 
P along (110) breaks up into multiple domains with local 
polarization along (100). We use these anomalies in the 
different energy contributions to identify Tq and Tg as 
function of strain and temperature. If no clear jump in 
energy is visible at Tc (due to the continuous character 
of this transition under strain) we instead use the 
lowest temperature with P = 0, which, however, gives 
rise to a large uncertainty of about 20 K. 

Fig. [3] (a) illustrates the so-obtained phase diagram of 
BTO under uniform (rj a = rjb) biaxial strain. In qualita¬ 
tive agreement with previous work , 13 ' 16 we find two tran¬ 
sition lines, Tc, corresponding to the appearance of in¬ 
plane and out-of-plane polarization, respectively, which 
cross for zero strain at the transition temperature of the 


























3 



UJ 

< 



T (K) 


FIG. 2. (Color online) (a) Domain structure under 0.75% tensile epitaxial strain for Tc > T = 400 K > Ts . Local dipole 
moments have been averaged along c and over 20 snapshots taken with a time distance of 0.2 ps. (b) Cross section of the 
polarization profile along [110]. Black: total polarization; Blue: Px; Red: P||, see text. Dashed line: fit to the domain profile 
given by Eq. ©■ (c) Temperature dependence of selected energy contributions for r/ = 0.75 %, normalized to 1 formula unit 
(f.u.) of BTO and taking the corresponding energy at 200 K as reference. Solid lines: case with multi-domain states (32x32x32 
cell); dashed lines: mono-domain case (16x16x16 cell). E\\ coupling between soft mode and local strain; Eu: local soft mode 
energy and inhomogeneous part of the elastic energy. Vertical lines indicate Tc and Ts. 


free bulk material. An additional transition line, Ts, 
separating multi-domain and single-domain phases is ob¬ 
served, extending from the tensile strain region into the 
region of small compressive strain. We note that Tg in 
Fig. E](a) is still somewhat dependent on the used cell size 
(e.g. for r\ = 0.75%, Ts is reduced by 35 K when using a 
92x92x92 cell, while the the domain wall width is fully 
converged for a 48x48x48 cell). 

Below both in-plane and out-of plane Tc, we obtain a 
multi-domain region where the in-plane components of 
the polarization exhibit the same domain patterns as in 
the tensile strain region (see e.g. Fig.[2](a)), but with an 
additional uniform polarization component along c, re¬ 
sulting in 60° domain walls parallel to {110}. Thus, in 
this region the local and global polarizations are along 
(110) and (111), respectively. We note that down to 
ri = —0.45%, the imposed strain acts as tensile strain 
relative to the shorter lattice constant of the tetrago¬ 
nal phase in the free material. For stronger compressive 
strain, Tg merges with the out-of-plane Tq and no multi- 
domain states are found. 

The case of general biaxial strain (with r] a ^ r]b) is il¬ 
lustrated in Fig. [3] (b)-(g). Under cooling, first a local 
polarization along ( 100 ) (light, red regions) appears in 
most regions of the phase diagram, with P pointing to¬ 
wards the longest lattice direction (a or b for r) a or 775 > 0 , 
c for compressive strain). The corresponding out-of-plane 
Tc decreases linearly with r] a +r]b, while the correspond¬ 
ing in-plane Tc increases linearly with the strain along 
the polarization direction, but is rather insensitive to the 
strain in the perpendicular direction. If two directions 
are under tensile strain, e.g. b and c for r) a = —0.85% 
and r]b = 0.35%, a multi-domain state with local polar¬ 
ization along (100) and 90° domain walls is more favor¬ 
able compared to a mono-domain state with polarization 
along (110). Thus, both mono- or multi-domain phases 


with local polarization along ( 100 ) are found below Tc- 

Under further cooling, local polarization along (110) 
(darker, blue) and then along ( 111 ) (black) becomes more 
favorable, analogously to the case of the free material. In 
general, the tendency for polarization along a certain di¬ 
rection increases, if the corresponding lattice constants 
in the clamped case are larger compared to their cubic, 
tetragonal, or orthorhombic counterparts in the free ma¬ 
terial. Thus, the transition into the (110) phases first sets 
in when both rj a and r)b are positive and for a combina¬ 
tion of large compressive and weak tensile in-plane strains 
(leading to an elongation along c). Again, multi-domain 
states occur for a broad range of strains and tempera¬ 
tures. For local (110) polarization, 90°, 60°, and 120° do¬ 
main walls are in principle possible^ As discussed above, 
only 60° domain walls parallel to { 110 } with local ac and 
be polarization are observed for r) a = rjb- For r] a ^ r/b, we 
find also 90° domain walls parallel to {100} with local 
ab/ab or ab/ab polarization and vanishing average global 
polarization along the shorter clamped lattice constant. 
For local (111) polarization, no multi-domain states are 
found for r) a = rjb- For r) a ^ rjb, 109° domain walls paral¬ 
lel to (010) and 71° walls parallel to (101) are possible^ 
and both types are found in our calculations. 71° walls 
occur for a combination of compressive and tensile strain 
and the corresponding local abc/abc or abc/abc polar¬ 
ization results in a vanishing global polarization along 
the compressed lattice constant. 109° walls with local 
abc/abc or abc/abc polarization, with neither net polar¬ 
ization along the shorter clamped lattice constant nor 
along c, are found for non-uniform tensile in-plane strain. 

In summary, we have shown, using ab initio -based 
molecular dynamics simulations, that the phase diagram 
of BTO under biaxial strain shows a variety of multi- 
domain regions. The formation of domains allows the 
system to recover the same sequence of para- and ferro- 
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FIG. 3. (Color online) Phase diagram of BaTiC >3 under biaxial strain. Colors/grayscale encode the paraelectric (white) and 
ferroelectric phases with local polarization along (100) (red/light gray), (110) (blue/dark gray), and (111) (black). Multi- 
domain configurations are indicated by striped regions as specified in the legend, (a) “Uniform” biaxial strain (r/a = ??(,). Open 
(filled) squares indicate out-of-plane (in-plane) Tc, triangles indicate Ts . The horizontal dashed line indicates Tc for the free 
bulk material, (b)-(g) General biaxial strain with r] a ^ r/b for different temperatures. 


electric phases as in the unstrained bulk material (at least 
for most combinations of rj a and 774 ,), while simultaneously 
fulfilling all epitaxial constraints. Such multi-domain 
states have not been found in previous ab initio -based 
studies of strained BTO, due to restrictions in the size 
of the simulation cells. However, similar multi-domain 
configurations, have been observed in previous empirical 
phase-field simulations^ Our results thus consolidate to 
some extent the empirical with the ab imfio-based sim¬ 
ulations, and provide new insights for future studies and 
the better interpretation of experimental data in strained 
BTO films. 
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